/* 
 * File:   main.cpp
 * Author: adriano
 *
 * Created on October 3, 2012, 5:49 PM
 */

#include <cstdlib>
#include <math.h>
#include <iostream>
#include "XmlFileReader.h"
#include "Entities/SimulationSetup.h"
#include "MeepExecutor.h"
#include <meep.hpp>
#include <vtk-5.10/vtkSmartPointer.h>
#include <vtk-5.10/vtkPolyDataMapper.h>
#include <vtk-5.10/vtkActor.h>
#include <vtk-5.10/vtkRenderer.h>
#include <vtk-5.10/vtkRenderWindow.h>
#include <vtk-5.10/vtkRenderWindowInteractor.h>
#include <vtk-5.10/vtkProperty.h>
#include <vtk-5.10/vtkPolyData.h>
#include <vtk-5.10/vtkCellLocator.h>
#include <vtk-5.10/vtkTriangleFilter.h>
#include <vtk-5.10/vtkXMLPolyDataReader.h>
#include <vtk-5.10/vtkButterflySubdivisionFilter.h>
#include <vtk-5.10/vtkLinearSubdivisionFilter.h>
#include <vtk-5.10/vtkUnstructuredGrid.h>
#include <vtk-5.10/vtkDataArray.h>
#include <vtk-5.10/vtkCellData.h>
#include <vtk-5.10/vtkAppendFilter.h>
#include <vtk-5.10/vtkDataSetMapper.h>
#include <vtk-5.10/vtkDataSet.h>
#include <vtk-5.10/vtkPoints.h>
#include <vtk-5.10/vtkIdList.h>
#include <vtk-5.10/vtkModifiedBSPTree.h>
#include <vtk-5.10/vtkOBBTree.h>
#include <vtk-5.10/vtkModifiedBSPTree.h>
#include <vtk-5.10/vtkGenericCell.h>
#include <vtk-5.10/vtkDecimatePro.h>

using namespace std;
using namespace meep;

/*
 * 
 */

vtkSmartPointer<vtkCellLocator> locator2D;
vtkSmartPointer<vtkOBBTree> locator3D;
vtkSmartPointer<vtkPolyData> polyData;
vtkSmartPointer<vtkIdList> idListOfIntersectionsInY;
vtkSmartPointer<vtkPoints> intersectionPointsInY;
double* objectBounds;
double TOLERANCE;

// variável que determina o centro do círculo
const double circleCenter[3] = {2., 0., 0.};
const double radius = 2.;

void print(SimulationSetup* setup)
{ 
    cout << "Dimensão: " << setup->GetDimension() << endl;
    cout << "\n";
    // computacional domain properties
    cout << " --Dominio Computacional -- " << endl;
    cout << "coordenada x inicial: " << setup->GetDomain()->GetCoordXInitial() << endl;
    cout << "coordenada x final: " << setup->GetDomain()->GetCoordXFinal() << endl;
    cout << "coordenada y inicial: " << setup->GetDomain()->GetCoordYInitial() << endl;
    cout << "coordenada x final: " << setup->GetDomain()->GetCoordYFinal() << endl;
    cout << "coordenada z inicial: " << setup->GetDomain()->GetCoordZInitial() << endl;
    cout << "coordenada z final: " << setup->GetDomain()->GetCoordZFinal() << endl;
    cout << "extensao em x: " << setup->GetDomain()->GetXExtension() << endl;
    cout << "extensao em y: " << setup->GetDomain()->GetYExtension() << endl;
    cout << "extensao em z: " << setup->GetDomain()->GetZExtension() << endl;
    cout << "resolucao do dominio: " << setup->GetDomain()->GetResolution() << endl;
    cout << "largura da PML: " << setup->GetDomain()->GetPmlThickness() << endl;
    cout << "\n";
    
    // source properties
    cout << "-- Fonte Eletromagnetica -- " << endl;
    cout << "Tipo de fonte: " << setup->GetSource()->GetType() << endl;
    cout << "Frequencia: " << setup->GetSource()->GetFrequency() << endl;
    cout << " Largura da frequencia: " << setup->GetSource()->GetWidth() << endl;
    cout << "Tempo inicial: " << setup->GetSource()->GetStartTime() << endl;
    cout << "Tempo final: " << setup->GetSource()->GetEndTime() << endl;
     cout << "\n";    
    
    // field properties
    cout << " --Campo -- " << endl;
    cout << "Tipo de campo: " << setup->GetField()->GetFieldType() << endl;
    cout << "Direção do campo: " << setup->GetField()->GetDirection() << endl;
    cout << "Posição inicial em x: " << setup->GetField()->GetInitialPosition()[0] << endl;
    cout << "Posicao final em x: " << setup->GetField()->GetFinalPosition()[0] << endl;
    cout << "Posição inicial em y: " << setup->GetField()->GetInitialPosition()[1] << endl;
    cout << "Posicao final em y: " << setup->GetField()->GetFinalPosition()[1] << endl;
    cout << "Posição inicial em z: " << setup->GetField()->GetInitialPosition()[2] << endl;
    cout << "Posicao final em z: " << setup->GetField()->GetFinalPosition()[2] << endl;
     cout << "\n";
}

void print(ProjectSettings* projectSettings)
{
    cout << "\n";
    cout << " -- Configuraçãoes do projeto -- " << endl;
    cout << "Diretório do projeto corrente: " << projectSettings->GetProjectDirectory() << endl;
    cout << "Projeto inclui filme?: " << projectSettings->IsMovie2D() << endl;
    cout << "\n";
}

vtkSmartPointer<vtkPolyData> readMaterialByXml(std::string file)
{    
    vtkSmartPointer<vtkXMLPolyDataReader> pdReader = vtkSmartPointer<vtkXMLPolyDataReader>::New();
    pdReader->SetFileName(file.c_str());
    pdReader->Update();
    
    return pdReader->GetOutput();
}
 
vtkSmartPointer<vtkPolyData> Triangulate(std::string file)
{      
    vtkSmartPointer<vtkPolyData> pd = readMaterialByXml(file.c_str()); 

    vtkSmartPointer<vtkTriangleFilter> triangulator = vtkSmartPointer<vtkTriangleFilter>::New();
    triangulator->SetInput(pd);
    triangulator->Update();
    
    vtkSmartPointer<vtkLinearSubdivisionFilter> subdivisor = vtkSmartPointer<vtkLinearSubdivisionFilter>::New();
    subdivisor->SetNumberOfSubdivisions(2);
    subdivisor->SetInput(triangulator->GetOutput());
    subdivisor->Update();
    
    return subdivisor->GetOutput();
}

void Set2DLocator()
{
    locator2D = vtkSmartPointer<vtkCellLocator>::New();
    locator2D->SetDataSet(polyData);
    locator2D->BuildLocator();
}

void Set3DLocator()
{
    locator3D = vtkSmartPointer<vtkOBBTree>::New();
    locator3D->SetDataSet(polyData);
    //locator3D->SetTolerance(TOLERANCE);
    locator3D->BuildLocator();
}

void SetObjectBounds()
{
    objectBounds = polyData->GetBounds();
}

void SetIntersectionPointsInY()
{
    intersectionPointsInY = vtkSmartPointer<vtkPoints>::New();
}

void SetIdListOfIntersectionsInY()
{
    idListOfIntersectionsInY = vtkSmartPointer<vtkIdList>::New();
}

void BuildWaveGuide(std::string path)
{
    std::string fileName = path + "/sembrPolyData.vtp";
    polyData = Triangulate(fileName);
    SetObjectBounds();
}

bool isOutside2DBounds(const vec &p)
{
    if ((p.x() < (objectBounds[0])) || (p.x() > (objectBounds[1])) ||
        (p.y() < (objectBounds[2])) || (p.y() > (objectBounds[3]))) 
    {
         return true;
    }
  
    return false;
}

bool isOutside3DBounds(const vec &p)
{
    if ((p.x() < (objectBounds[0])) || (p.x() > (objectBounds[1])) ||
        (p.y() < (objectBounds[2])) || (p.y() > (objectBounds[3])) ||
        (p.z() < (objectBounds[4])) || (p.z() > (objectBounds[5])))
    {
         return true;
    }
  
    return false;
}

double epsilon2D(const vec &p)
{
    //cout << "coord x : " << p.x() << " - coord y: " << p.y() << "coord z " << p.z() << "\n";
    // outise bound object
    if (isOutside2DBounds(p)) return 1.0;
      
    // on surface object
    // it's seems that dimensions not set need to have its value equals zero
    double currentVec[3] = {p.x(), p.y(), 0.};   
           
    vtkIdType cellId = locator2D->FindCell(currentVec);
    //cout << "coord x : " << currentVec[0] << " - coord y: " << currentVec[1] << "\n";
    if (cellId >= 0  )
    {
//        if (polyData->GetCellData()->GetScalars()->GetTuple1(cellId) > 1.21)
//                cout << "Cell Id: " << cellId << " - Dieletric value : " << polyData->GetCellData()->GetScalars()->GetTuple1(cellId) << endl;        
        
        return polyData->GetCellData()->GetScalars()->GetTuple1(cellId);
    }
    
    return 1.0;
}

double epsilon3D(const vec &p)
{
    // outise bound object
    // return dieletric of air
    if (isOutside3DBounds(p)) return 1.0;     
    // initial and final points to build the ray intersection
    // test: using tolerance to avoid the initial point to be over waveguide bound
    const double initialYPoint[] = {p.x(), objectBounds[2] - TOLERANCE, p.z()}; 
    const double finalPoint[] = {p.x(), p.y(), p.z()};
    // locator doing intersection
    locator3D->IntersectWithLine(initialYPoint, finalPoint, intersectionPointsInY, idListOfIntersectionsInY);    
    // if there is some intersecion, the number of pointsInY will be different of zero
    // if the number of points is odd, the position vector is inside the waveguide
    if ((intersectionPointsInY->GetNumberOfPoints() != 0) && (intersectionPointsInY->GetNumberOfPoints() % 2 != 0))
    {
        int nrOfIds = idListOfIntersectionsInY->GetNumberOfIds();
        // get last intersection through ids intersection to use scalar value based on cell with that id
        vtkIdType idType = idListOfIntersectionsInY->GetId(nrOfIds - 1);
        //cout << "Dieletric Value: " << polyData->GetCellData()->GetScalars()->GetTuple1(idType) << endl;
        return polyData->GetCellData()->GetScalars()->GetTuple1(idType);
    }
    // return dieletric of air if terms above aren't reached
    return 1.0;
}

// teste 1: dielétrico do silício em um retângulo infinito X 1
// referência: http://ab-initio.mit.edu/wiki/index.php/Meep_Tutorial
double epsRectangule(const vec &p)
{
    // retângulo centrado na origem do espaço cartesiano - ponto(0, 0);
    // a altura tem limite inferior em -0.5 e limite superior em 0.5, ou seja, igual a 1
    //  não há limite na largura, esta é infinita
    if (p.y() >= -0.5 && p.y() <= 0.5)
    {
        return 12.;
    }
    
    return 1.;
}

// teste 2: dielétrico do sílicio em um círculo sólido
double epsCircle(const vec &p)
{
    double currentVec[3] = {p.x(), p.y(), 0.};
    
    double value = ((currentVec[0] - circleCenter[0]) * (currentVec[0] - circleCenter[0])) + 
                   ((currentVec[1] - circleCenter[1]) * (currentVec[1] - circleCenter[1]));
    
    if (value <= (radius * radius)) // está dentro ou na fontreira do círculo
        return 12.0; // valor dielétrico do silício
    
    return 1.0; // valor dielétrico do ar
}

// teste 3: dielétrico do vidro (silica) em uma esfera sólida
double epsSphere(const vec &p)
{   
    double value = ((p.x() - circleCenter[0]) * (p.x() - circleCenter[0])) + 
                   ((p.y() - circleCenter[1]) * (p.y() - circleCenter[1])) +
                   ((p.z() - circleCenter[2]) * (p.z() - circleCenter[2]));
    
    if (value <= (radius * radius)) // está dentro ou na fontreira da esfera
        return 2.25; // valor dielétrico do vidro
    
    return 1.0; // valor dielétrico do ar   
}

// teste para gerar arquivo para tese
double epsVoid(const vec &p)
{       
    return 1.0; // valor dielétrico do ar   
}

double epsFiberMonomode(const vec &p)
{
    // validate core region
    if ((p.x() >= -4.5 && p.x() <= 4.5) && (p.y() >= -0.41 && p.y() <= 0.41))
    {
        return 2.15561124;
    } // validate cladding region
    else if ((p.x() >= -4.5 && p.x() <= 4.5) && (p.y() >= -6.25 && p.y() <= 6.25))
    {
        return 2.14007641;
    }
    else // return air
    {
        return 1.0;
    }
    // infinite length in x direction
    // validate core region
//    if (p.y() >= -0.41 && p.y() <= 0.41)
//    {
//        return 2.15561124;
//    } // validate cladding region
//    else if (p.y() >= -6.25 && p.y() <= 6.25)
//    {
//        return 2.14007641;
//    }
//    else // return air
//    {
//        return 1.0;
//    }
}

void ExecuteSimulationInVoid(int argc, char** argv)
{
    meep::initialize mpi(argc, argv);
    
    std::string path = "/home/adriano/Projects/SEM-BR/trunk/sembr.meep/sembr_output/3DSimulation";
    const double pml_thickness = 0.1;
    const double x_extension = 1. + pml_thickness;
    const double y_extension = 1. + pml_thickness;
    const double z_extension = 1. + pml_thickness;
    const double resolution = 10;
    
    grid_volume grid_vol = vol3d(x_extension, y_extension, z_extension, resolution);
    grid_vol.center_origin();
    
    structure strct(grid_vol, epsVoid, pml(pml_thickness));
    strct.set_output_directory(path.c_str());
    
    continuous_src_time src(0.15, 6.67, 0, 10); 
    
    fields field(&strct);
    field.set_output_directory(path.c_str());
    field.use_real_fields();       
    field.output_hdf5(Dielectric, grid_vol.surroundings());    
    field.add_point_source(Ez, src, vec(-1, 0, 0));
        
    // outup
    h5file *h5_file = field.open_h5file("ez", h5file::WRITE, "", true);
        
    while (field.time() < field.last_source_time())
    {   
        field.output_hdf5(Ez, grid_vol.surroundings(), h5_file, true, false, "");
        field.step();        
    }
}

void ExecuteSimulationInRectangule(int argc, char** argv)
{
    meep::initialize mpi(argc, argv);
    
    std::string path = "/home/adriano/Projects/SEM-BR/trunk/sembr.meep/sembr_output/pureMeep_Rectangule/ez";
    const double pml_thickness = 1.;
    const double x_extension = 14. + pml_thickness;
    const double y_extension = 8. + pml_thickness;
    const double resolution = 10;
    
    grid_volume grid_vol = vol2d(x_extension, y_extension, resolution);
    grid_vol.center_origin();
    
    structure strct(grid_vol, epsRectangule, pml(pml_thickness));
    strct.set_output_directory(path.c_str());
    
    continuous_src_time src(0.15, 6.67, 0, 50); 
    
    fields field(&strct);
    field.set_output_directory(path.c_str());
    field.use_real_fields();       
    field.output_hdf5(Dielectric, grid_vol.surroundings());    
    field.add_point_source(Ez, src, vec(-7, 0));
        
    // outup
    h5file *h5_file = field.open_h5file("ez", h5file::WRITE, "", true);
        
    while (field.time() < field.last_source_time())
    {   
        field.output_hdf5(Ez, grid_vol.surroundings(), h5_file, true, false, "");
        field.step();        
    }
}

void ExecuteSimulationInCircle(int argc, char** argv)
{
    meep::initialize mpi(argc, argv);
    
    std::string path = "/home/adriano/Projects/SEM-BR/trunk/sembr.meep/sembr_output/pureMeep_Circle";
    const double pml_thickness = 1.;
    const double x_extension = 14. + pml_thickness;
    const double y_extension = 6. + pml_thickness;
    const double resolution = 30;
    
    grid_volume grid_vol = vol2d(x_extension, y_extension, resolution);
    grid_vol.center_origin();
    
    structure strct(grid_vol, epsCircle, pml(pml_thickness));
    strct.set_output_directory(path.c_str());
    
    continuous_src_time src(0.25, 4., 0, 100); 
    
    fields field(&strct);
    field.set_output_directory(path.c_str());
    field.use_real_fields();       
    field.output_hdf5(Dielectric, grid_vol.surroundings());
    
    field.add_point_source(Ez, src, vec(-2, 0));
    
    // outup
    h5file *h5_file = field.open_h5file("ez", h5file::WRITE, "", true);
    
    while (field.time() < field.last_source_time())
    {   
        field.output_hdf5(Ez, grid_vol.surroundings(), h5_file, true, false, "");
        field.step();        
    }
}

void ExecuteSimulationInSphere(int argc, char** argv)
{
    meep::initialize mpi(argc, argv);
    
    std::string path = "/home/adriano/Projects/SEM-BR/trunk/sembr.meep/sembr_output/pureMeep_Sphere";
    const double pml_thickness = 0.5;
    const double x_extension = 6. + pml_thickness;
    const double y_extension = 4. + pml_thickness;
    const double z_extension = 3. + pml_thickness;
    const double resolution = 30;
    
    grid_volume grid_vol = vol3d(x_extension, y_extension, z_extension, resolution);
    grid_vol.center_origin();
    
    structure strct(grid_vol, epsSphere, pml(pml_thickness));
    strct.set_output_directory(path.c_str());
    
    continuous_src_time src(0.15, 6.67, 0, 100); 
        
    fields field(&strct);
    field.set_output_directory(path.c_str());
    field.use_real_fields();       
    field.output_hdf5(Dielectric, grid_vol.surroundings());    
    field.add_point_source(Ex, src, vec(-2., 0., 0.));
        
    // outup
    h5file *h5_file = field.open_h5file("ex", h5file::WRITE, "", true);
        
    while (field.time() < field.last_source_time())
    {   
        field.output_hdf5(Ex, grid_vol.surroundings(), h5_file, true, false, "");
        field.step();        
    }   
}

void ExecuteSimulationInFibermonomode(int argc, char** argv)
{
    meep::initialize mpi(argc, argv);
    
    std::string path = "/home/adriano/Projects/SEM-BR/trunk/sembr.meep/sembr_output/fiber_monomode";
    const double pml_thickness = 1.;
    const double x_extension = 18. + pml_thickness;
    const double y_extension = 18. + pml_thickness;
    const double resolution = 10;
    
    grid_volume grid_vol = vol2d(x_extension, y_extension, resolution);
    grid_vol.center_origin();
    
    structure strct(grid_vol, epsFiberMonomode, pml(pml_thickness));
    strct.set_output_directory(path.c_str());
    // 193.423THz
    gaussian_src_time src(193.423 * pow(10, 11), 0.155); 
    
    fields field(&strct);
    field.set_output_directory(path.c_str());
    field.use_real_fields();       
    field.output_hdf5(Dielectric, grid_vol.surroundings());    
    //field.add_point_source(Ez, src, vec(-2.5, 0));
    field.add_volume_source(Ey, src, volume(vec(-4.5, -0.5), vec(-4.5, 0.5)));
    field.add_volume_source(Hz, src, volume(vec(-4.5, -0.5), vec(-4.5, 0.5)));
        
    // outup
    h5file *h5_file = field.open_h5file("ey", h5file::WRITE, "", true);
        
    while (field.time() < field.last_source_time())
    {   
        field.output_hdf5(Ey, grid_vol.surroundings(), h5_file, true, false, "");
        field.step();        
    }
}

void ExecuteSimulationInMicrolens(int argc, char** argv) 
{
    meep::initialize mpi(argc, argv);
    string executablePath = argv[0];
    int execPathLength = executablePath.length();
    executablePath = executablePath.substr(0, (execPathLength - 36));
    
    cout << endl;
    cout << "executable path: " << executablePath << endl; 
    string projectPath = executablePath + "/project_settings/";
    // read setup from xml file
    XmlFileReader xmlReader;
    // project settings config
    ProjectSettings* psPtr = xmlReader.getProjectSettings(projectPath + "project_settings.xml"); // TODO: temporary for testing. Need to be changed!
    print(psPtr);
    // simulation config
    SimulationSetup* setupPtr = xmlReader.GetSetup(psPtr->GetProjectDirectory() + "/pre_processing/setup.xml");
    print(setupPtr);    
    // construct waveguide based on vtkPolyData defined on simulation
    BuildWaveGuide(psPtr->GetProjectDirectory() + "/pre_processing");
    Set3DLocator();
    SetIntersectionPointsInY(); 
    SetIdListOfIntersectionsInY();
    std::string post_proc_path = psPtr->GetProjectDirectory() + "/post_processing";
    // setup meep
    grid_volume grid_vol = vol3d(setupPtr->GetDomain()->GetXExtension(), 
                                 setupPtr->GetDomain()->GetYExtension(), 
                                 setupPtr->GetDomain()->GetZExtension(),
                                 10);
    grid_vol.center_origin();
    
    structure strct(grid_vol, epsilon3D, pml(setupPtr->GetDomain()->GetPmlThickness()));
    strct.set_output_directory(post_proc_path.c_str());
    // 193.423THz
    gaussian_src_time src(193.423 * pow(10, 11), 0.155); 
    
    fields field(&strct);
    field.set_output_directory(post_proc_path.c_str());
    field.use_real_fields();  
    field.output_hdf5(Dielectric, grid_vol.surroundins()); 
    field.add_volume_source(Ey, src, volume(vec(-0.5, 6.74, -0.5), vec(0.5, 6.74, 0.5)));
    //field.add_volume_source(Hz, src, volume(vec(-0.5, 6.74, -0.5), vec(0.5, 6.74, 0.5)));
        
    // outup
    h5file *h5_file = field.open_h5file("ey", h5file::WRITE, "", true);
        
    while (field.time() < field.last_source_time())
    {   
        field.step();        
    }
    
    field.output_hdf5(Ey, grid_vol.surroundings(), h5_file, true, false, "");
}

void ExecuteVisualization(int argc, char** argv)
{
    string executablePath = argv[0];
    int execPathLength = executablePath.length();
    executablePath = executablePath.substr(0, (execPathLength - 36));
    
    cout << endl;
    cout << "executable path: " << executablePath << endl; 
    string projectPath = executablePath + "/project_settings/";
    // read setup from xml file
    XmlFileReader xmlReader;
    // project settings config
    ProjectSettings* psPtr = xmlReader.getProjectSettings(projectPath + "project_settings.xml"); // TODO: temporary for testing. Need to be changed!
    print(psPtr);
    // simulation config
    SimulationSetup* setupPtr = xmlReader.GetSetup(psPtr->GetProjectDirectory() + "/pre_processing/setup.xml");
    print(setupPtr);    
    // construct waveguide based on vtkPolyData defined on simulation
    BuildWaveGuide(psPtr->GetProjectDirectory() + "/pre_processing");
    
    vtkSmartPointer<vtkDataSetMapper> ds = vtkSmartPointer<vtkDataSetMapper>::New();
    
    if (locator2D != NULL) {
        Set2DLocator();
        ds->SetInput(locator2D->GetDataSet());
    }        
    else {
        Set3DLocator();
        ds->SetInput(locator3D->GetDataSet());
    }
                
    ds->SetScalarModeToUseCellData();
       
    vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
    actor->SetMapper(ds);
    //actor->GetProperty()->SetColor(0., 0., 0.);
    
    vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
    renderer->AddActor(actor);
    double backgroundColor[3] = {255., 255., 255.};    
    renderer->SetBackground(backgroundColor);
    
    vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
    renderWindow->AddRenderer(renderer);
    
    vtkSmartPointer<vtkRenderWindowInteractor> interactor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
    interactor->SetRenderWindow(renderWindow);
    
    renderWindow->Render();
    interactor->Start();
}

void ExecuteMeep(int argc, char** argv)
{
    // work with path to sembr.meep executable
    // taking off the name of executable file in argv[0]
    // the information needed here is only the path to sembr.meep project
    string executablePath = argv[0];
    //string systemPath = argv[1];
    int execPathLength = executablePath.length();
    // TODO: trocar 36, abaixo, por 11
    executablePath = executablePath.substr(0, (execPathLength - 36));
   
    cout << endl;
    //cout << "system path: " << systemPath << endl;
    cout << "executable path: " << executablePath << endl; 
    //string projectPath = systemPath + executablePath + "/project_settings/";
    string projectPath = executablePath + "/project_settings/";
          
    // read setup from xml file
    XmlFileReader xmlReader;
    // project settings config
    ProjectSettings* psPtr = xmlReader.getProjectSettings(projectPath + "project_settings.xml");
    //print(psPtr);
    // simulation config
    SimulationSetup* setupPtr = xmlReader.GetSetup(psPtr->GetProjectDirectory() + "/pre_processing/setup.xml");
    //print(setupPtr);   
    
    TOLERANCE = 1 / setupPtr->GetDomain()->GetResolution();
    // construct waveguide based on vtkPolyData defined on simulation
    BuildWaveGuide(psPtr->GetProjectDirectory() + "/pre_processing");
    MeepExecutor* meepExec;   
    
    // 3D simulation
    if (setupPtr->GetDimension() == 3)
    {
        Set3DLocator();
        SetIntersectionPointsInY(); 
        SetIdListOfIntersectionsInY();
        meepExec = new MeepExecutor(argc, argv, psPtr, setupPtr, epsilon3D);
        meepExec->Run();
    }
    else // 2D simulation
    {        
        Set2DLocator();
        meepExec = new MeepExecutor(argc, argv, psPtr, setupPtr, epsilon2D);
        meepExec->Run();
    }
    
    cout << "\n\n";
    cout << "MPI elapsed time: " << meepExec->GetElapsed_time() << endl;
}

int main(int argc, char** argv) 
{   
    //ExecuteSimulationInCircle(argc, argv);
    //ExecuteSimulationInRectangule(argc, argv);
    //ExecuteSimulationInVoid(argc, argv);
    //ExecuteSimulationInSphere(argc, argv);
    //ExecuteSimulationInFibermonomode(argc, argv);
    ExecuteSimulationInMicrolens(argc, argv);
    //ExecuteVisualization(argc, argv);    
    //ExecuteMeep(argc, argv);

    return 0;
}
